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ABSTRACT 

The astronomical theory of the ice ages is investigated using a simple climate model which 
includes the ice sheets explicitly. A one-level, zonally averaged, seasonal energy-balance equation 
is solved numerically for sea-level temperature T" as a function of latitude and month (similar to 
North, 1975). Seasonally varying snow cover (which affects planetary albedo) is included 
diagnostically by parameterizing monthly snowfall and snowmelt in simple ways. The net annual 
accumulation and ablation on the ice sheet surface at each latitude are computed using the same 
parameterizations as for snow cover above (with T corrected for ice sheet height using a lapse 
rate of — 6.5°C km"'). Treatment of the ice sheets follows Weertman (1976) with ice flow 
approximated as perfect plasticity, which constrains the ice sheet profiles to be parabolic. The 
northern hemisphere’s ice sheet is constrained to extend equatorward from 75 ®N (corresponding 
to the Arctic Ocean shoreline). 

Model ice age curves are generated for the last several 1(X) Kyears by computing the seasonal 
climate as above once every 2 Kyears, with insolation calculated from actual Earth orbit 
perturbations. The change in ice sheet size for each 2 Kyear time step depends only on the net 
annual, snow budget integrated over the whole ice sheet surface. In these model runs, the 
equatorward tip of the northern hemisphere’s ice sheet oscillates through ~7° in latitude, 
correctly simulating the phases and approximate amplitude of the higher frequency components 
(-'43 Kyear and 22 Kyear) of the deep-sea core data (Hays et al., 1976). However, the model 
fails to simulate the dominant glacial-interglacial cycles (~100 to 120 Kyear) of this data. The 
sensitivity of the model ice age curves to various parameter changes is described, but none of 
these changes significantly improve the fit of the model ice age curves to the data. In the 
concluding section we generalize about the types of mechanisms that might yield realistic 
glacial-interglacial cycles. 


1. Introduction and summary 

There is little agreement as yet on the dominant 
causes of the Quaternary ice ages, although many 
mechanisms have been suggested (described in 
Beckinsale, 1973; Andrews, 1975, p. 71). Data 
from deep-sea sediment cores provide continuous 
records of some climatic variables over the last 
several 100 Kyears (e.g., Broecker and Van Donk, 
1970; Hays et al., 1976; Emiliani, 1978). These 
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records of global ice sheet volume have been 
interpreted as showing quasi-periodic glacial-inter- 
glacial cycles with fast retreats from maximum to 
minimum volumes occurring at intervals of approx- 
imately 100 to 120 Kyears; superimposed on these 
cycles are secondary oscillations with smaller 
amplitudes and higher frequencies. Others have 
cautioned that some or all of the fluctuations in 
these records may not be periodic but essentially 
random (e.g., Shackleton, 1969, p. 145; Lemke, 
1977). In any case, these continuous records are 
suitable for comparisons against any simulated 
records generated by quantitative models for- 
mulated to test the various suggested ice age 
mechanisms. 
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One such mechanism, the “astronomical” or 
“Milankovitch” theory, involves variations of 
incoming solar insolation due to the secular 
perturbations of the Earth’s orbit. These variations 
are the only well-known external forcing of climate 
on ice age time scales. The historical development 
of the astronomical theory is lucidly described in 
Imbrie and Imbrie (1979). Milankovitch estimated 
the sensitivity of climate to this forcing using 
insolation-curve computations (e.g., Milankovitch, 
1941), and these insolation-curves have been 
refined and extended by Van Woerkom (1953), 
Vernekar (1972), Berger (1978) and others. Other 
recent investigations have used seasonal climate 
models (Shaw and Donn, 1968; Budyko and 
Vasishcheva, 1971; Saltzman and Vernekar, 1971; 
Suarez and Held, 1976, 1979; Schneider and 
Thompson, 1979), explicit ice sheet models 
(Weertman, 1976; Birchfield and Weertman, 

1978) , a combined seasonal cUmate-ice sheet 
model (Pollard, 1978), and a combined annual 
mean climate-ice sheet-ocean model (Sergin, 

1979) . Also Calder (1974) and Imbrie and Imbrie 
(1980) have used single non-explicit “response 
equations”. Most of these studies have generated 
simulated ice age curves of one sort or another, and 
comparisons with deep-sea core records con- 
sistently indicate that the astronomical forcing can 
account for the observed secondary oscillations; 
this positive result is consistent with the power 
spectrum analyses of Hays et al. (1976) and 
Kominz and Pisias (1979) (but see Evans and 
Freeland, 1977). However, none of the models 
above have correctly simulated the dominant 
glacial-interglacial cycles of the records as a 
response to the astronomical forcing. 

This paper reports on a combined seasonal 
climate-ice sheet model that was described briefly 
in a preliminary paper (Pollard, 1978). We had 
hoped that the additional non-linearity due to the 
interactions between the ice sheets and the seasonal 
cycle might be the missing factor required to 
produce realistic glacial-interglacial cycles. With 
this motivation we explored the model’s sensitivity 
to a systematic range of parameter variations. In 
Section 2 below, the model is formulated and its 
solution for the present seasonal climate is 
described. In Section 3, the model’s long-term 
response to the orbital perturbations is presented, 
and in Section 4 we describe the sensitivity of this 
response to small changes in various parameter 


values and types of parameterizations. In the; 
sections we emphasize some differences betwec 
these sensitivities and those of some other clima; 
models, and suggest which differences in tf 
models can account for the different sensitivitie 
Although our parameter changes have relative! 
slight effects on the fit to the present climate, son 
of them have significant effects on the mode! 
simulated ice age curves, i.e., on its long-ten 
response to the astronomical forcing. Over most c 
the parameter ranges these curves retain th 
secondary oscillations observed in the deep-se 
core records; however, with the model in its presen 
form we are still unable to generate any curve 
resembling the observed dominant glacial-inter 
glacial cycles. 

Therefore, this paper supports the results of thi 
earlier studies mentioned above. Further, since thi 
present model still cannot account for the glacial- 
interglacial cycles, we suggest in Section 5 that ont 
or more other long-term mechanisms may b« 
important. 


2. Model formulation and present-day 
results 


The model has previously been outlined in 
Pollard (1978). There are two distinct parts corres- 
ponding to the two distinct time-scales of the global 
seasonal climate and the long-term ice sheet 
response. In Sections 2.1, 2.2 and 2.3 the seasonal 
part of the model and its solution for the present 
climate are described, and in Section 2.4 the ice 
sheets are incorporated into the model. 


2.1. Seasonal energy-balance equation 
Following North (1975), the climate through one 
year over a spherical globe is described by a 
zonally averaged, one-level energy-balance 
equation for sea-level air temperature T-. 




D 8 


+ [A-t-BT\ = Q{xm-a) + S (1) 


Here x is sin (latitude) and t is time. All dependent 
variables are defined as ~1 month running-means, 
so daily correlations are effectively assumed cons- 
tant. Boundary conditions are (1 — dTidx = 
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0 at the poles x = ±1 (Hantel, 1972). Other 
symbols and their values for our “standard” model 
are: 

C = 4.6 X 10’ J m~’ °C“', a constant seasonal 
heat capacity of the atmosphere-land-ocean 
system (equivalent to a layer of liquid water 
1 1 m thick). 

D = 0.501 ,x 10‘ m’ s“', a linear diffusion 
coefficient acting over the whole thickness of 
the layer represented by C. 

R = 6.36 X lO* m, radius of the Earth. 

A = 207 W m"’, B = 1.9 W m“’ °C~‘, net infrared 
radiation coefficients. 

Q = zonal mean insolation at the top of the 
atmosphere, computed for any given era from 
the orbital elements in Berger (1976) using a 
solar constant of 1360 W m~’. 
a = ra^ + (1 — r) Of, earth-atmosphere albedo, 

= 0.62 represents areas covered by seasonal 
snow or ice sheet (see below), and Of = 0.3 1 + 
0.08 l(3x’ — 1)/21 represents areas free of 
seasonal snow and ice sheet. We set r = 1 
north of 75 °N to represent perennial Arctic 
Ocean sea-ice (Ku and Broecker, 1967; 
Hunkins et al., 1971), and r = 1 south of 70 ° S 
to represent a fixed Antarctic ice sheet (see 
below). At all other latitudes r = 0.6 when 
covered by seasonal snow or ice sheet, and r = 
0 when free of seasonal snow and ice sheet. 

S = 1.27 (J m“’ s“'l per (g cm"^ month”*], 
representing latent heat of fusion released or 
required at each latitude by the varying 
amount of seasonal snow cover (see below). 
The annual mean of S at each latitude is zero. 
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uncertainty exists for the prediction of cloud 
amounts even in much more complex models. 

2.2. Seasonal snowmelt and snowfall 

Seasonally varying snow cover on land at 
sea-level is modelled diagnostically by parameter- 
izing monthly snowmelt and snowfall as functions 
of the current air temperature T and insolation Q. 
We use 

Snowmelt (g cm"’ month"') = max 10; aT{°C) 

+ i)Q(W m-’) -p cl (2) 

Equation (2) is basically an energy-balance 
equation for a melting snow/ice surface with 
seasonal heat storage neglected, and as such is 
equivalent to the snowmelt parameterization of 
Suarez (1976). Equation (2) is also used below for 
the monthly ablation on ice sheet surfaces (Section 
2.4). For the standard model we use a = 10, b = 
0.32, c = —47; these values and the adequacy of 
this parameterization for the ice age problem are 
discussed in more detail in Pollard (1980). 

Whereas snowmelt is a micrometeorological 
process, snowfall depends on synoptic-scale pro- 
cesses and the basic dependence of past and 
seasonal variations of snowfall on the zonally 
averaged variables T and Q is not nearly so 
apparent as for snowmelt. For this reason [and 
following Suarez (1976)] we set 

Snowfall (g cm"’ month"') 

fPdat.) ifr<0°C 

= (3) 

[ 0 ifr>0°C 


Most of these parameterizations [discussed for 
instance in Coakley (1979)] have found general use 
in many annual mean energy-balance models and 
are based on annual mean data. There is currently 
some doubt about the physical basis of the diffusive 
heat parameterization (Van Loon, 1979), but it 
may be valid for global and seasonal scales 
(Hartjenstein and Egger, 1979; Lorenz, 1979). The 
infrared radiation and albedo parameterizations are 
found to represent seasonal data less accurately 
than annual mean data, possibly due to independ- 
ently varying cloud cover (White, 1976; Warren 
and Schneider, 1979), but we use them here for 
simplicity. We also neglect possible variations of 
cloud cover in past eras, which might be serious for 
the ice age problem. Unfortunately, considerable 


where P is the present observed zonal and annual 
mean precipitation rate at each latitude. The 
northern hemispheric data given in Schutz and 
Gates (1971—4) is used for F(lat.) in both model 
hemispheres, since the present precipitation in 
south polar regions is clearly affected by Antarctic 
topography. Equation (3) is also used for the 
monthly accumulation on ice sheet surfaces in the 
standard model (Section 2.4), but some effects of 
ice sheet topography on the local accumulation rate 
will be included later. 

For simplicity the model neglects variations of 
sea-ice. Although this might be serious for the ice 
age problem, seasonal and past variations of sea-ice 
in the northern hemisphere are somewhat smaller 
(by factors of ~ 2 to 3) than those of seasonal snow 
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cover and ice sheets on land (Flint, 1971, figs. 4.8, 
4.9; Emiliani and Getss, 1957, table 2; Saltzman 
and Vernekar, 1975, table 1). 

2.3. Fit to present climate 

As in North (1975) and North and Coakley 
(1979), the solution of eq. (1) is simplified by 
expressing the latitudinal dependences of T and the 
right-hand forcing in terms of the eigenfunctions of 
the spherical diffusion operator. Only the first three 
eigenfunctions are kept [i.e., Legendre polynomials 
1, X, and (3 a:^ — l)/2], since any finer latitudinal 
resolution would probably not be realistic due to 
the coarse parameterizations used. However, the 
diagnostic snow budget parameterizations (2) and 
(3) prevent a corresponding convenient Fourier 
expansion in time t (cf. North and Coakley, 1979). 
Equation (1) is numerically integrated forward in 
time-steps of 1 month through consecutive years 
until initial transients decay to negligible levels and 
a repeated seasonal cycle is attained (usually after 
- 10 years). 

Fig. 1 shows the sea-level temperature solution 
of the standard model, for the present orbital 
elements and with no northern hemispheric ice 
sheet. For comparison table la in Warren and 
Schneider (1979) shows the equivalent data for 
surface air temperature. The parameter values of 
the standard model given above were chosen to fit 
the present temperature data of the northern 
hemisphere. With the infrared radiation, albedo and 
snow budget parameterizations all fixed (some- 
what arbitrarily), C was adjusted to yield realistic 
seasonal amplitudes at mid and high latitudes and 
D was adjusted to yield realistic latitudinal 
gradients. 

The most obvious discrepancy in Fig. 1 from 
reality is the excessive seasonal amplitudes in the 
southern hemisphere due to the uniform value of C; 
in this respect our climate model is effectively two 
northern hemispheres patched together. The 
seasonal snow cover in the model southern hemi- 
sphere is only a crude analogy for the seasonal 
variation of sea-ice around Antarctica. Perhaps 
more seriously for the application to the Milan- 
kovitch theory, the model seasonal cycle in both 
hemispheres is lagged ~2 months behind the 
insolation cycle, which is about ~1 month more 
than observed in the northern hemisphere. This 
excessive lag has also been found by North and 
Coakley (1979) and Thompson and Schneider 
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Fi%. I . Zonal mean sea-level temperature vs latitude and 
month for standard model present-day solution. Values 
shown are degrees centigrade. The dashed curve shows 
the latitudes of maximum temperature at each month, 
with maximum temperature values for some months in 
parentheses. The dotted curves show the latitudinal 
extents of seasonal snow cover in each hemisphere, for 
the months when snow exists equatorward of the limits of 
the Arctic Ocean (75 °N) or Antarctic ice sheet (70“S). 


(1979) in their simplest model versions, and is 
related to the lack of any longitudinal contrast in 
seasonal heat capacity between land and ocean. 
This excessive lag is reflected in the seasonal snow 
line shown in Fig. 1; the onset of snow in autumn is 
~1 month later than observed (Kukla, 1975). 

Fig. 2 plots the net annual “potential” snow 
budget for the northern hemisphere corresponding 
to the present-day solution. This shows the net 
annual snowfall minus snowmelt that would occur 
on any mountain glacier or ice cap surface at a 
given latitude and elevation h, calculated for each 
month from (2) and (3) but with T corrected to T— 
6.5(°C km~')-^ to allow for the atmospheric lapse 
rate. (Of course, snowmelt can continue below the 
nominal level for these surfaces so that negative net 
annual budgets are possible.) The zero-budget line 
generally agrees with present data such as the 
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Fig. 2. Net annual snow accumulation minus ablation 
that would occur in the model present-day solution on 
hypothetical ice surfaces at given latitudes and 
elevations. Values shown are g cm"^ year~‘. A model ice 
sheet profile given by (4) roughly representing Greenland 
is also shown, although its contribution to the albedo a in 
(1) was ignored for the present-day solution. The dashed 
line is the observed “regional snowline” of Paschinger 
(1912), redrawn from Sugden and John (1975, Fig. 5.5). 


generalized “regional snowline” (Paschinger, 
1912), glacier altitudes in western U.S.A. (Meier, 
i960), and the altitude of the equilibrium line in 
southern Greenland (Schuster, 1954). However, the 
model zero-budget line is ~500 m too high in polar 
latitudes north of ~70° N. There are few estimates 
of the present net annual budget for the entire 
Greenland ice sheet; approximately half of the total 
mass loss is by the calving of icebergs (Paterson, 
1969, p. 228). The net annual accumulation minus 
ablation for the model Greenland profile in Fig. 2 
averages to —4 g cm“^ over the southern half and 
-f21 g cm“^ over the northern half. 

2.4. Ice sheets 

Ice sheets are incorporated into the climate 
model following Weertman’s (1964, 1976) simple 
treatment. Ice sheet flow under its own weight is 
approximated to be perfectly plastic, which con- 
strains the model ice sheet profiles to always 
remain parabolic; 


h(s) = U(L-lsl)]‘'^ 


(4: 
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where h is the elevation of the ice sheet surfaci 
above sea-level, L is the half-width and s is thi 
latitudinal distance from the ice sheet center (with i 
taken positive towards the equator). The ice sheet 
base is assumed to be isostatically depressed tc 
depths 0.5 A(s) below sea-level. 2 is a constant 
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proportional to the yield stress of ice; for our 
standard model we use A = 10 m, corresponding to 
a yield stress of ~0.7 bars. This value is slightly less 
than in Weertman (1976) but still gives somewhat 
greater central thicknesses (by ~30%) than those 
modelled by Paterson (1972) and those suggested 
by Greenland and Antarctica today. 

One model ice sheet, representing the Laurentide 
and Scandinavian ice sheets of past eras, is 
constrained to extend equatorward with its 
northern tip fixed at 75“ N (corresponding to the 
Arctic Ocean shoreline). Where the margins of the 
real northern hemispheric ice sheets reached 
continental coastlines, further advance was pre- 
vented by rapid iceshelf and iceberg calving into the 
ocean, but their equatorial extent and overall 
volume were probably limited more by ablation on 
their southern flanks (cf. Flint, 1971, p. 484, 600). 
Therefore, as in Weertman (1976), the long-term 
variation in the model ice sheet size is controlled by 
the net accumulation (snowfall) minus ablation 
(mostly snowmelt) on its southern half only (i.e., s 
> 0). Also, since its profile is constrained by (4), 
any change in size is determined simply by the total 
ice volume added to or removed from the entire 
southern half. Writing the net annual accumula- 
tion minus ablation (per unit surface area) as m{s), 
then the change in ice sheet size per unit long-term 
time increment dr is given by: 


dV 

dr 


= — 

dz 


1 mjs) 
pJo (1 year) 


(5) 


V is the volume of the southern half of the ice sheet 
(per unit longitudinal distance), related to L from 
(4). p is the mean ice sheet density, taken as 0.9 g 
cm“’. 

The other model ice sheet representing An- 
tarctica is circular and centered on the South Pole. 
The existing Antarctic ice sheet is prevented from 
advancing beyond the continental bedrock by rapid 
calving into the southern oceans, so we do not 
allow this model ice sheet to ever extend equator- 
ward beyond 70° S. In fact, the ice sheet remains at 
this maximum size in all ice age runs shown below, 
since its net annual budget is always positive. 
Significant ice age fluctuations in the real Antarctic 
ice sheet might have occurred (e.g., Wilson, 1964), 
but we defer this possibility to future model 
developments. 

Model ice age curves are generated by solving 
(1) for one year’s “climate” T(x,t) once at the start 
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of each long-term time-step dr, using the current 
orbital elements. The presence of the ice sheets 
influences T only through the albedo a. The 
monthly snowmelt and snowfall at any point on the 
ice sheet surface is given by (2) and (3) with T 
corrected to T— 6.5(°C and with h given 

by (4). These values are integrated to give the net 
annual snow budget for the northern hemispheric 
ice sheet which is then used in (5) to determine the 
change in its size over the next dr years. For all ice 
age curves below, we used dr = 2 Kyears. This 
type of procedure involving two distinct climatic 
time scales is described more formally in Hassel- 
mann (1976), and has been applied to the ice age 
problem by Eriksson (1968) using a generalized 
analytical approach. 

3. Ice age results: Standard model 

3.1. Ice sheet response 

The ice age simulation produced for the last 400 
Kyears by our standard model is compared in Fig. 
3 with a (5 ’*0 deep-sea core record from Hays et al. 
(1976). The model curve of northern hemispheric 
ice sheet volume reflects the forcing of both 
obliquity (~41 Kyear period) and precession (~22 
Kyear period), with obliquity dominating in times 


of small eccentricity (e.g., 60 to 0 Kyears BP), 
response of an ice sheet model to indivi 
sinusoidal forcings with these periods has 1 
described by Birchfield (1977).] The model 
sheet curve closely resembles an equivalent 
sheet curve in Weertman (1976, his fig. 6) and 
the response-curves in Imbrie and Imbrie (19' 
and coincides both in phase and approxim 
amplitude with the secondary oscillations of 
<5'®0 record. The long-term volume inertia of l 
ice sheets has produced a phase lag of ~5 to 
Kyears behind the orbital forcing (by vis 
comparison with various insolation-curves), co 
sistent with the measured phase relationships 
Hays et al. (1976). But despite its explicit season 
climate treatment, our model still lacks the dom 
nant glacial-interglacial cycles and the drastic ic 
sheet retreats (e.g., 18 to 0 Kyears BP), in commo 
with Weertman’s and Birchfield’s models. 

In all ice age runs the pattern of net annua 
accumulation minus ablation over the ice sheets 
remains the same as that shown for Greenland in 
Fig. 2, i.e., widespread net accumulation over the 
central regions and much stronger net ablation over 
a relatively narrow strip near the equatorward tip. 
Therefore, the non-linear possibility sugeested by 
Weertman (1964), of retreating ice sheets becoming 



Fig. 3. (a) 5 “O recorded from two combined deep-sea cores (from ~45° S lat.), redrawn from Hays et al. (1976, Fig. 
9). (b) Ice age curve for standard model, showing northern hemispheric ice sheet volume normalized by its volume 
with equatorward tip at 50” N. Right-hand scale shows corresponding latitudes of its equatorward tip. The dashed 
curves from 100 to 0 Kyears BP show the effect of choosing different initial ice sheet sizes at 100 Kyears BP. (c) 
Maximum (monthly mean) sea-level temperature at 55°N (solid curve) and 55°S (dashed curve), corresponding to the 
ice age run in (b). 
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acnant and thus having a shorter shrinkage 
mc-scale than growth time-scale, does not occur 
anv of our model runs. Imbrie and Imbrie’s 
1980) "response-equation” does contain such a 
on linearity and produces some suggestions of 
lacial-inierglacial cycles and drastic ice sheet 
ctreais in response to the eccentricity variations. 
Consequently their fit to the ^ '*0 record is slightly 
better than ours, but stills leaves room for 
improvement. 

The effect of orbital changes on the standard 
model can be viewed in another way in Fig. 4. Each 
curve shows how the northern hemispheric ice 
sheet would grow or recede as a function of its size, 
for various combinations of the orbital elements. 
This figure is a development of fig. 7 in Weertman 
(1961), who analyzed the general shape of the 
curves as a simple geometrical consequence of ice 
sheet topography and the typical net snow budget 
pattern. As expected (Berger, 1978), the grouping 
of the curves in Fig. 4 shows that eccentricity plus 
precession are most important for large (low 
latitude) ice sheets, whereas obliquity has an 
important effect for small (high latitude) ice sheets. 
For fixed orbital elements, the ice sheet size would 
move along the relevant curve until it either arrives 
at a stable point S or is ablated away to the Arctic 
shoreline. The ice sheet tip in the ice age run of Fig. 


3 (b) varies between ~50°N and ~60“N, and this 
is the region in Fig. 4 where approximately half of 
the orbital combinations are enlarging the ice sheet 
towards stable points 5 and the other half are 
ablating it toward the Arctic shoreline. However, 
this potential amplitude of some 30° in latitude 
implied by Fig. 4 is reduced to ~10° latitude in Fig. 
3 (b) by the volume inertia of the ice sheet. 

The full amplitude of ~30° in latitude is reflected 
in the dashed curves in Fig. 3 (b), which show the 
effect of choosing different initial ice sheet sizes at 
100 Kyears BP. If the initial tip position is north of 
~61°N, the ice sheet does not return to the solid 
curve. Instead it is ablated back to the Arctic 
shoreline and remains there forever, since the 
majority of orbital combinations produce negative 
regimes for these small ice sheet sizes. This is an 
example of an “intransitive” climate system (Lorenz, 
1970), with two possible stable branches depend- 
ing on the initial conditions. The orbital forcing by 
itself is insufficient to produce transitions from one 
branch to another, and so we effectively had to 
“choose” the ice-sheet-free branch for the present- 
day fit in Section 2.3, and the other branch 
exhibiting ice sheet oscillations for this section. 
Some modifications described below will make 
these transitions possible (Fig. 10) or will diminate 
the distinction between the two branches (Fig. 8). 
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Fig. 4. Curves of net annual northern hemispheric ice sheet budget (averaged over southern half) as functions of its 
size, for various combinations of the orbital elements. (Precession is defined here so that it is zero when the northern 
hemispheric winter solstice coincides with perihelion, and 180° when this solstice coincides with aphelion.) These 
curves are for the standard model. For fixed orbital elements, “S” are stable equilibrium points, “U" are unstable 
equilibrium points, and the arrows show directions of ice sheet growth and decay. 
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The time-scale on which the dashed curves in 
Fig. 3 (b) converge onto one branch or another is 
on the order of '•40 Kyears. This represents the 
intrinsic time-scale of the quasi-linear ice sheet 
response on either branch, and is roughly con- 
sistent with the amount oflag (mentioned above) of 
the solid curve in Fig. 3 (b) behind the orbital 
forcing. Because the model has only one long-term 
time-derivative (in (5)1, there can be no free internal 
oscillations of the system in the absence of external 
forcing, as there are in Sergin’s (1979) climate- 
ocean-ice sheet model. 


3.2. Temperature response 

Fig. 3 (c) shows the variation of maximum 
summer temperature at particular latitudes for the 
standard ice age run. The temperature curve for 
55° N seems (visually) to lag only ~0 to 4 Kyears 
behind the orbital forcing, considerably less than 
the ice sheet volume curve. This summer temper- 
ature curve has roughly the same phase as various 
northern hemispheric curves in Shaw and Donn 
(1968), Suarez and Held (1979) and Schneider and 
Thompson (1979). The higher-frequency compo- 
nent of the temperature curve for 55° S reflects the 
forcing of precession, which is 180° out of phase 
between the two hemispheres. The phase of the 
55° S curve does not agree with that of the 
secondary oscillations of the “Ts” deep-sea core 
record for ~45° S in Hays et al. (1976); the latter is 
found to lead those of the ^‘*0 core record by ~2 


Kyears. This disagreement may be due to tli 
shortcomings in the model formulation for th 
southern hemisphere. [Both Suarez and Held'- 
(1979) and Schneider and Thompson’s (1979 
models have realistic southern hemispheric 
seasonal heat capacities, but they appear to yield 
southern hemispheric temperature phases that are 
opposite to each other.] 

Table 1 summarizes the sensitivity of northern 
hemispheric temperatures to some particular 
changes in the orbit and ice sheet size (analogous to 
Fig. 4). The most extreme orbital variation with no 
northern hemispheric ice sheet (col. 4, Table 1) 
causes seasonal temperature changes comparable 
to those caused through the albedo a by a full 
glacial-interglacial ice sheet variation at fixed orbit 
(col. 5, Table 1). However this ice sheet variation 
causes annual mean temperature changes at mid 
and high latitudes that are considerably larger than 
those caused by the orbital variations, in agree- 
ment with the trend in Sergin (1979, fig. 16). 

We now compare the sensitivities in Table 1 to 
those of other models. The sensitivities of annual 
mean temperatures to the orbital variations (cols. 2 
to 4, Table 1) are similar to those of Schneider and 
Thompson’s (1979) model, but are generally less 
than 1/4 of those in Suarez and Held (1979). The 
extra sensitivity in the latter model is probably due 
to greater albedo feedback of seasonal snow. As 
discussed by Suarez and Held, albedo feedback is 
greater for models like theirs having realistic 
land-ocean longitudinal contrast. However in 


Table 1. Differences of sea-level temperatures at particular latitudes between different orbits {as defined 
in Fig. 4) with no northern hemispheric ice sheet {columns 2 to 4), and between different ice sheet sizes 
with the same orbit {column 5) 


Latitude 


6 orbit d orbit d orbit 



d ice sheet tip 
(73°N-50°N) 


89° N 


1.1 

(-ii) 

-0.4 

U’f) 

0.7 


45° N 

u.:) 

0.1 

ao 

-0.3 

U:l) 

-0.1 


1°N 


-0.7 

U?) 

0.0 

(4:9 

-0.7 



The values in parentheses are differences in °C for the months of maximum northern hemispheric summer 
temperatures (upper value) and minimum northern hemispheric winter temperatures (lower value). The other value 
outside the parentheses is the annual mean temperature difference in °C. 
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addition to this effect, Suarez and Held have 
adjusted their parameter values so that over nearly 
the whole range of orbital variations, their northern 
hemispheric minimum summer snowline remains 
equatorward of the Arctic Ocean; in contrast, our 
model solutions (for cols. 2 to 4, Table 1) have no 
seasonal snow (i.e., no albedo feedback) in summer 
equatorward of our Arctic Ocean shoreline at 
75° N, as in Fig. 1. The real summer snowline 
seems to be roughly intermediate between these two 
model situations (Dickson and Posey, 1967; 
Kukla, 1975; WilUams, 1978). 

The reductions of T in summer due to the ice 
sheets (col. 5, Table 1) are comparable to those 
found by CLIMAP (1976) data for sea-surface 
temperatures of the last glacial maximum. The 
corresponding air temperature reductions found by 
the more complex zonal model of Saltzman and 
Vernekar (1975) and the GCM of Gates (1976a,b) 
are generally larger than in col. 5 (by ~ 100%), 
whereas the GCM of Williams (1974) has found 
much larger air temperature reductions (of 
;^20°C). 

These differences in the temperature sensitivities 
of the various models lead to some uncertainty in 
applications to the Milankovitch theory. For 
instance, it may be that if we replaced our simple 
climate model with Suarez and Held’s model, the 
amplitude of the ice sheet response in Fig. 3 (b) (as 
well as the annual mean temperature response) 
would be increased by a factor of ~4. However, 
our simple model already realistically simulates the 
secondary oscillations of the data, and any 
large alteration of its ice sheet sensitivity would 
destroy this agreement. What is needed is a model 
modification to produce the dominant glacial- 
interglacial cycles of the data. 

4. Ice age results: parameter sensitivity 

The aim of this section is to explore the range of 
parameter values and parameterizations that still 
yield realistic ice age secondary oscillations, and 
also to search for a model modification that may 
produce the full glacial-interglacial cycles. We do 
not concentrate on the sensitivity of particular 
climate solutions T(x,t) per se. Most of the 
parameter changes examined below do slightly 
perturb the present-day solution of Section 2.3 
(e.g., by ^±2°C in Fig. 1, by ^±500 m in the 
zero-budget line altitude in Fig. 2), but these 


perturbations are minor compared to the existing 
coarseness of the fit to present data. 

4.1. Climate parameter variations 

Fig. 5 shows the sensitivity of the ice age curve 
to small variations in the diffusion parameter D, 
with all other parameters held constant. Only the 
response for the last 100 Kyears is shown, but this 
is sufficient to demonstrate the behavior for any 
time period. We have performed the same exercise 
for small variations (-2 to 10%) in each one of the 
“climate parameters” C, B, Of, a^, a, b, and c 
[appearing in (1) and (2)] in turn, all with virtually 
the same result as in Fig. 5. 

For parameter changes that have the effect of 
reducing the summer temperatures in mid and high 
latitudes or reducing ablation for a given temper- 
ature and insolation, the mean position of the ice 
age curves move smoothly equatorward. The effect 
of precession (^22 Kyear period) becomes more 
pronounced equatorward of ~45°N, as might be 
expected from Fig. 4. Apart from this effect, the 
basic model response is unchanged and still 
resembles the secondary oscillations of the ice age 
data. [For larger parameter variations than shown 
here, a sudden transition to ice-sheet-covered 
northern hemispheric continents might be expected. 
In other runs (not shown) investigating solar 
constant variations, ice sheets could exist in stable 
equilibrium down to ~20° N for solar constant 
reductions of up to ~5%, beyond which they grew 
down to the equator.] For parameter changes in the 
other direction, there is a sudden transition to an 
ice-sheet-free northern hemisphere; there is no 
stable mean ice sheet tip position between •~55°N 
and the Arctic shoreline. Once the minimum size is 
reached, the ice sheet never grows out again [as in 
Fig. 3 (b)]. 

The behavior in Fig. 5 can be explained by 
referring to Fig. 4. Climate parameter variations 
basically have the effect of shifting the pattern of 
curves in Fig. 4 horizontally relative to the 
“regime” x-axis. For small shifts toward more 
negative regimes, ice sheet sizes in the range ~50° 
to 60° N begin to lie more on the unstable ablating 
branches of most orbital curves, and so are reduced 
back to the Arctic shoreline. For relatively large 
shifts toward more positive regimes, more orbital 
curves have stable points S and the range of 
latitudes of these stable points smoothly moves 
equatorward. 
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Fig. S. Ice age curves for the standard model except: 
(a) Diffusion coefficient D = (0.31/0.30) x standard 
model value, (b) Diffusion coefficient D = (1) x standard 
model value, (c) Diffusion coefficient D = (0.29/0.30) x 
standard model value, (d) Diffusion coefficient D — 
(0.28/0.30) X standard model value, (e) Diffusion 
coefficient D = (0.23/0.30) x standard model value. 


If only one climate parameter is varied at a time, 
the range of interest is limited by the sudden retreat 
of the ice sheets las in Fig. 5 (a)). Another common 
sensitivity test (Coakley, 1979; Warren and 
Schneider, 1979) is to vary two or more para- 
meters simultaneously so that their basic effects 
partially cancel each other and the mean ice sheet 
position can be held “on scale”. Three examples of 
this are shown in Fig. 6. In curves (a) and (b) large 
variations in D are compensated by changes in the 
albedo contrast to maintain nearly the same 
temperature field for a given orbit and ice sheet 
size. In curves (c) a low value of the seasonal heat 
capacity C (which produces larger seasonal cycles 
and higher summer temperatures at high latitudes) 
is compensated crudely by reducing the ablation 
rate by a constant factor. 



Fig. 6. Ice age curves for the standard model except; 
(a) Solid curve: Diffusion coefficient D = (0.23/0.30) x 
standard model value, and % = o)** + 0.017, = 

0 */' — 0.14, where and o*'* are the standard model 
albedos for snow-ice/free and snow-ice/covered surfaces 
respectively, (b) Dashed curves: Diffusion coefficient D 
= (0.37/0.30) X standard model value, and = a)*' 
— 0.015, Bj = a*** + 0.134, showing two different choices 
of initial ice sheet size, (c) Dotted curves: Seasonal 
heat capacity C = (7/11) x standard model value, and 
with ablation rates (2) reduced by a factor 0.42, 
showing two different choices of initial ice sheet size. 

For the large value of D and the low value of C, 
two different initial ice sheet sizes are chosen to 
show that the lower branch of stable ice sheet tip 
positions has shifted south to ~40° to 45° N 
latitude. This is because both these parameter 
variations favor the net snow budget of large ice 
sheets relative to small ones, due to increased ice 
sheet albedo feedback for large B and due to larger 
summer temperature increases at higher latitudes 
for low C. Consequently the stable points “S” in 
the corresponding orbital-curve diagrams (not 
shown) are located ~10° to 15° further south than 
in Fig. 4. This is the only real difference in Fig. 6 
from the standard model response, and the ampli- 
tude and phase of the ice sheet oscillations have 
remained basically the same. 

We have not repeated the exercise in Fig. 6 for 
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Z all possible combinations of parameter variations 
§ lof which there are on the order of 10’, as noted by 
I Imbrie and Imbrie (1980)1. In the development of 
I this model we experimented relatively unsystem- 
{ atically with many (perhaps ~200) combinations of 
I parameter values that appear in eqs. (1) to (5), and 
I never found any types of ice age response other 
^ than those described in this paper. 

1 The constants in the ablation parameterization 
(2) are not tightly constrained by present glacial 

r data (Pollard, 1980). We ran several ice age curves 
f (not shown) using widely different values; for 

2 instance (a,b,c) = (20,0,80) respectively so that 
ablation depended only on temperature, and (a,b,c) 

^ = (5,0.32,-68) so that the temperature depen- 
dence was half that of the standard model. In these 
2 runs, the responses to the orbital perturbations 
were basically unchanged from the standard model, 
producing secondary oscillations of the same phase 
and magnitude without any suggestion of full 
glacial-interglacial cycles. [The temperature depen- 

- dence in (2) cannot be eliminated completely. With 
a = 0, the model yields unrealistic seasonal cycles 

i with snow-free high latitudes in summer and 

■ perennial snow in mid-latitudes, due to the 
latitudinal forms of the precipitation rate and the 
seasonal insolation forcing. In the present model, 
ablation and seasonal snow cover must be control- 
led mainly by temperature.) 

Henderson-Sellers and Meadows (1979) and 
Cogley (1979) have suggested that variations of 
high-latitude ice cover have affected the planetary 
j albedo to a much lesser degree than in many other 
models. Correspondingly Fig. 7 (a) shows a run 
. with no albedo feedback at all, neither from the 
seasonal snow nor from the ice sheets; a is simply a 
t: constant function of latitude. The present-day 
solution with this parameterization still fits the 
present data as well as the standard model in 

- Section 2.3. The ice age response in Fig. 7 (a) is 
; basically unchanged from that of the standard 
, model, suggesting that ice-albedo feedback is not a 
^ significant mechanism for the secondary oscil- 

■ lations of the ice age records. Fig. 7 (b) shows 
* another run using an albedo partly dependent on 
1' solar zenith angle, as investigated by Lian and Cess 

(1977). Again there was no basic change in the ice 
age response, as might be expected from the 
indifference of Fig. 7 (a) to the albedo para- 
meterization. 

h Fig. 7 (c) shows an ice age run for an annual 



Fig. 7. Ice age curves for the standard model except: (a) 
Solid curve: Using fixed albedos, a = 0.35 + 0.2[(3x’ — 
l)/2l always, (b) Dashed curve: Using fraction of hourly 
insolation absorbed by the earth-atmosphere system = 
/jO.35 cos (s) -t 0.479], where s is hourly solar zenith 
angle and / = 1 or 0.6 for snow/ice-free or snow/ice- 
covered surfaces respectively, (c) Dotted curve: Using 
annual mean insolation in eq. (1), and using (o,6,c) = 
(10,0, 122) respectively in eq. (2). Also setting S = 0 in 
eq. (I). 

mean version of our model. For this version, the 
current annual mean insolation at each latitude is 
used for Q in (1), and ablation depends only on T. 
The resulting annual mean climate solutions for 
individual years are much the same as in North 
(1975), with a sea-level snowline potentially at 
~-70°N (except that for Fig. 7 (b), this latitude 
region is occupied by ice sheet). The peak-to-peak 
amplitude of the ice sheet tip response of the annual 
mean version is reduced to ~1.5° in latitude, and 
temperature variations at fixed latitudes are all 
^1 °C. Therefore we find that the seasonal cycle is 
necessary for our model to produce realistic 
secondary oscillations. We have seen that seasonal 
albedo feedback is not necessary for the seasonal 
model’s ice age response [Fig. 7 (a)], so seasonal 
albedo feedback cannot be the important difference 
between the seasonal and annual mean models. The 
important difference seems to be due to the fact 
that the orbital perturbations change the seasonal 
cycles of temperature (at the latitudes around the 
ice sheet tip) much more than the annual mean 
temperatures. The ice sheets respond just as much 
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to changes in the seasonal cycles as to the annual 
mean changes, due mostly to the non-linearity in 
the ablation parameterization (2), and so the 
seasonal version of the model produces a much 
larger ice sheet response. 

We now compare the sensitivities in Fig. 7 to 
those of two other seasonal energy-balance models. 
North and Coakley’s (1979) model has a seasonal 
snowline, a perennial “ice line” fixed to the — 10°C 
zonal and annual mean isotherm, and also con- 
tains longitudinal land-ocean asymmetry. In res- 
ponse to obliquity changes of 1.2°, their perennial 
ice line changes by 3° in latitude for the seasonal 
version and 2° in latitude for the annual mean 
version. Presumably for obliquity variations of 2.4° 
(more representative of the actual orbital pertur- 
bations) this would imply ice line changes of 6° lat. 
(seasonal version) and 4° lat. (annual mean 
version). In contrast, our model ice sheet tip varies 
by 7° lat. (seasonal version) and ~1.5° lat. (Mnual 
mean version) in response to the actual orbital 
perturbations (including precession). This points to 
an important difference between the two models: 
their ice line responds only to the mean annual 
temperature and so the increased response of their 
seasonal version is due to slight variations in the 
residual correlation between the seasonal cycles of 
albedo and insolation. However, as discussed 
above, our more non-linear ice sheets respond 
directly to variations in the seasonal cycles them- 
selves, resulting in a greater difference in response 
between seasonal and annual mean versions. 

Schneider and Thompson (1979) find that the 
sensitivity of temperature to the orbital pertur- 
bations in their seasonal climate model is decreased 
by '•30% by using fixed (constant) albedos 
compared to using seasonally varying albedos. 
When our model is run with no ice sheets we find 
basically the same result as theirs, and the different 
result implied by Fig. 7 (a) is due to the presence of 
the ice sheets. The albedos in our standard model 
can change only in the winter months when the 
seasonal snowline extends equatorward beyond the 
ice sheet tip; during the summer months the ice 
sheets prevent any change in albedo and so the 
seasonal variation of albedo is reduced consider- 
ably (cf. discussion of Table 1). 

North and Coakley (1979) and Thompson and 
Schneider (1979) both find that the sensitivities of 
their models to 1% solar constant variations are 
nearly the same for seasonal and annual mean 


versions, which at first sight contradicts the r 
in Fig. 7 (c). However, solar constant variat 
primarily affect the annual mean insolation and 
the seasonal cycles, and so are a fundamem 
different type of forcing from the orbital per 
bations in Fig. 7. [In fact we do find (not sho 
that the sensitivities to small solar consi 
variations of our seasonal and aimual m 
versions are very nearly the same, with glo 
annual mean temperatures changing by 1.8 °C 
1 % change in solar constant.) 

We now mention two other modifications to 
model that were tried. In some runs, monthly zo 
precipitation was parameterized as a function 
sea-level temperature T and ^TV^latitude, 
opposed to the fixed precipitation of the standa 
model. Several similar functions were tried, f| 
instance 


Precip. (g cm“^ month'*) 

= max [4;l 10 dT{°C)!d lat. (deg.)l) 
exp[r(°C)/17] 

This is a very rough fit to present seasonal zon: 
data in Schutz and Gates (1971-4); simila 
parameterizations have been investigated b; 
Schneider and Thompson (1977). The functioi 
implies a reduction in precipitation during glacia 
maxima associated with lower saturation vapoi 
pressures, which has sometimes been suggested ai 
a significant ice age factor. However, these para 
meterizations produced no ice age curves signifi 
cantly different from those of the standard model 
suggesting that the secondary oscillations of the ic< 
sheet records have been caused more by ablatioi 
variations than accumulation variations. Data ii 
Yapp and Epstein (1979) and Ruddiman anc 
McIntyre (1979) are suggestive of important ic( 
age precipitation variations due to changini 
longitudinal land-ocean temperature contrasts, bu 
this is outside the scope of the present model. 

In several runs (not shown) we crudely at 
tempted to simulate the long-term effect of th( 
short-term “random” weather variability, as 
analyzed by Hasselmann (1976) and Lemk< 
(1977). In these runs a random term, rectangularlj 
distributed between ±10 g cm"’ year"* [i.e., ±2 x 
10^ g cm"’ (2 Kyear)"*], was added to the meat 
ice sheet budget at each 2 Kyear time step. Th( 
resulting ice sheet volume curves were not signifi 
cantly different from the standard model response 
with no suggestion of any drastic ice sheet retreats. 
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4.2. Ice sheet parameter variations 

No really different ice age responses in either the 
amplitude or the phase of the secondary oscil- 
lations have been produced by the climate para- 
meter variations above. As shown below, para- 
meter variations concerning the ice sheet can have 
somewhat greater effects. 

Fig. 8 show ice age runs with the precipitation on 
all ice sheet surfaces reduced by a factor 
exp l-/i(km)/3j from the zonal mean, where h is 
the local ice sheet elevation given by (4). To 
balance this, ablation is also reduced slightly. This 
crudely models the topographic blocking of storms 
carrying precipitation to the ice sheet interiors, as 
observed on Antarctica and Greenland today 
(Mock, 1967; Chorlton and Lister, 1968). The 
effect of slight variations in the weather para- 
meters in Fig. 8 is similar to Fig. 5, but now there is 
no sudden transition to an ice-sheet-free northern 
hemisphere, and stable mean ice sheet tip positions 
can exist between ~55° N and the Arctic shoreline. 
The corresponding orbital-curve diagram is shown 
in Fig. 9; these curves do not bend back to negative 
regimes for small ice sheets nearly so much as in 
Fig. 4, allowing small stable ice sheets. This also 
implies that ice age runs with different initial ice 
sheet sizes converge to the same curve, as shown 
by the dashed curve in Fig. 8; the equivalent curve 
in Fig. 3 (b) retreated to the Arctic shoreline. 



Fig. 9. As Fig. 4 except that precipitation on ice sheets is reduced by factor exp I— h(km)/3] from the zonal mean, and 
b = 0.26 in eq. (2). 



Fig. 8. Ice age curves for the standard model except that 
precipitation on ice sheets is reduced by factor 
exp 1— A(km)/3I from the zonal mean, and: (a) Ablation 
coefficient b = 0.30 in eq. (2). (b) Ablation coefficient 
b = 0.28 in eq. (2). (c) Ablation coefficient b = 0.26 in eq. 
(2). (d) Ablation coefficient b = 0.24 in eq. (2). (e) 
Dashed curve is as for curve (c) but with different initial 
ice sheet size, (f) Dotted curve is for standard model 
except that precipitation on ice sheets is altered by factor 
2 exp 1— h(km)/3] from the zonal mean. 

The change from Fig. 4 to Fig. 9 can be 
explained as follows: the reduction in precipitation 
is greater for large ice sheets than for small ones, 
whereas the ablation reduction we have used to 
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balance these reductions afTects all ice sheet sizes 
equally. Any equivalent modification to the model, 
that favors the net snow budget of small ice sheets 
relative to large ice sheets, produces changes in the 
same direction as Fig. 9; for instance, smaller lapse 
rate magnitudes than 16.51 °C km~’ in (2), thinner 
ice sheet profiles, or less ice sheet albedo feedback. 
(The trend in the opposite direction was described 
briefly in connection with Fig. 6.) 

In contrast to the precipitation reductions in ice 
sheet interiors, the steep flanks of ice sheets can 
locally increase precipitation on the sides facing the 
prevailing winds (Mock, 1967). To crudely test this 
effect we ran some ice age curves (not shown) with 
precipitation on all ice sheet surfaces increased by a 
factor of 2 over the zonal mean (and with ablation 
similarly increased by a constant factor). How- 
ever, the only effect on the response was the 
predictable one of doubting the peak-to-peak 
amplitude of the secondary oscillations to ~15° in 
latitude. There were still no suggestions of a 
realistic glacial-interglacial cycle. 

The amplitude of the response can also be 
increased by using thinner ice sheets, i.e. by 
reducing the value of A in eq. (3) below 10 meters. 
Fig. 10 shows two runs using A » 4 meters, which 
lies slightly below the range of values appropriate 
for existing ice bodies (Paterson, 1972, table 2). 
Perhaps more unrealistically, isostatic depression 
of the land surface beneath the ice sheet is ignored. 
These curves show transitions between a mean ice 
sheet position around ~55°N (from 300 to 200 
and from 100 to 0 Kyears BP), and a much smaller 
mean position trapped near the Arctic shoreline 


(from 2(X) to 100 Kyears BP). This type of 
response is intermediate between those in Figs. 5 
and 8; now the increased amplitude of the basic 
oscillations (forced by obliquity and precession) is 
sufficient to ocassionally bridge the gap between 
the two stable positions. This new situation might 
be classified as “almost intransitive” (Lorenz, 1970). 

We have chosen slightly different minimum ice 
sheet sizes for the two curves in Fig. 10, and this 
difference can occasionally be important in allow- 
ing an “escape” from the Arctic shoreline or not 
(e.g., at ~150 Kyears BP). In reality regional land 
topography becomes important for these nascent 
ice sheets (Loewe, 1971; Barry et al., 1975), and 
the present model just suggests that such details for 
small ice sheets can sometimes affect the form of 
the subsequent response. 

TTie curves in Fig. 10 are notably similar to 
many of the curves produced by the ice sheet 
models of Weertman (1976) and Birchfield and 
Weertman (1978), and also to the curve produced 
by Calder’s (1974) response equation. The ice sheet 
thicknesses used by Birchfield and Weertman (A = 

14 meters) are more similar to our standard model 
value, and considerably greater than those in Fig. 
10. Presumably this is compensated by our climate 
model producing smaller net snow budget varia- 
tions (due to the orbital forcing) than those 
produced by their more geometrical parameter- 
ization; the curve in Weertman (1976), Fig. 6) that 
is most similar to our standard model curve [Fig. 3 ' 
(b)l uses accumulation and ablation values that are 
generally ~l/2 to ~l/3 of those in our model. As | 
in their studies, we have also tried using circular 



Fig. 10. Ice age curves for standard model except with no isostatic depression below the ice sheets, and; (a) Solid 
curve: A = 4.2 meters in eq. (3), and values of (a,b,c) in eq. (2) are (8,0.234,-39.2) respectively. Also minimum ice 
sheet half-width = 1.5° lat. (b) Dashed curve: A = 3.6 meters in eq. (3), and values of (a,b,c) in eq. (2) are 
(8,0.226,-39.2) respectively. Minimum ice sheet half-width = 2.0° lat. 


Tellus 32 (1980), 4 







■ars BP). This type of 
oetween those in Figs. 5 
amplitude of the basic 
■:iquity and precession) is 
bridge the gap between 
This new situation might 
ransitive” (Lorenz, 1970). 
iy different minimum ice 
irves in Fig. 10, and this 
>y be important in allow- 
- Arctic shoreline or not 
)■ In reality regional land 
;Jortant for these nascent 

• Barry et al., 1975), and 
■gests that such details for 
etimes affect the form of 

} are notably similar to 

• iduced by the ice sheet 
976) and Birchfield and 
Iso to the curve produced 
ise equation. The ice sheet 
tfield and Weertman (A = 
lar to our standard model 
greater than those in Fig. 
•mpensated by our climate 
r net snow budget varia- 
ital forcing) than those 
s geometrical parameter- 
rtman (1976), Fig. 6) that 
ndard model curve [Fig., 3 
nd ablation values that are 
of those in our model. As 
s also tried using circular 

D 

D 

5 

2 

Q. 

5 

cr « 

<i> I 

5 


he ice sheets, and; (a) Solid 
jectively. Also minimum ice 
:s of (a,b,c) in eq. (2) are 


Tellus 32 (1980), 4 





' 1 I I I - . ' ■ vt - 



RESPONSE OF MODEL TO ORBITAL PERTURBATIONS DURING QUATERNARY ICE AGES 315 


rather than linear northern hemispheric ice sheets, 
and have taken the ice sheet budget over the whole 
surface and not just the southern half. However 
these modifications leave our ice age runs basically 
unchanged. Although their curves and those in Fig. 
10 might be considered suggestive of possible 
longer-period cycles and their spectra may contain 
some power at periods > 100 Kyears (Birchfield 
and Weertman, 1978), the models still fail to 
produce realistic glacial— interglacial cycles. In fact, 
the “secondary” oscillations in Fig. 10 are much 
larger than those observed in the records. 

5. Concluding section 

In the various ice age runs above, the ampli- 
tudes of the northern hemispheric ice sheet volume 
fluctuations are generally between 20% to 50% of 
the maximum glacial volume, corresponding to 
between and 10° in the latitude of its 
equatorward tip. These fluctuations generally agree 
both in phase and amplitude with the secondary 
oscillations of (J’*0 deep-sea core records, at least 
to within the small variations from record to record 
and within the mixed effects of ice sheet volume 
and ocean temperature in the <5**0 signal (Emiliani 
and Shackleton, 1974). 

The components of the model that are necessary 
to produce this response are the explicit treatment 
of ice sheet topography and snow budget, and the 
seasonal and latitudinal variations of temperature. 
In fact, given these features we cannot find any 
reasonable parameter changes that do not give 
realistic secondary oscillations (except for cases 
giving an ice-sheet-free northern hemisphere). 
Using annual mean insolation in (1) reduces the ice 
age response by a factor of ~4, but eliminating 
albedo feedback has very little effect on the model’s 
response (see Fig. 7). However, albedo feedback of 
the ice sheets might still be important for full 
glacial-interglacial cycles. These sensitivities are 
related to the ice sheet and ablation parameteri- 
zations, as discussed in connection with Fig. 7. 

We have not been able to produce realistic 
glacial-interglacial cycles with the present model. 
Starting at small ice sheet size, the model can 
plausibly simulate the relatively slow ~80 Kyear 
growths to glacial maximum [for instance by 
adjusting the ice sheet precipitation parameteri- 
zation as in Fig. 8 (f)]; it is the drastic ~20 Kyear 


retreats back to interglacials that the model lacks. 
The estimates of Laurentide ice sheet volumes in 
Paterson (1972) imply mean ice sheet budgets 
averaging — 50 g cm"^ yr~* between 14 and 9 
Kyears BP. (Ice sheet retreat after this point was 
probably accelerated strongly due to being split by 
marine waters of Hudson Bay at ~8 Kyears BP.) 
In our model ice age runs the mean northern ice 
sheet budget varies only between ~ -t-20 and — 20 g 
cm'^ yr~', and this range is fairly independent of 
model ice sheet details (for instance, the full ice 
sheet retreat from 20 to 0 Kyears BP is achieved 
artificially in Fig. 10 by reducing the ice sheet 
volume inertia, but the mean ice sheet budget in this 
period is still — 20 g cm"^ yr“'). What can change 
in the ice sheet environment to produce mean 
budgets of around —50 g cni~^ yr"' between 14 and 
9 Kyears BP, and also produce generally positive 
budgets ybr e/te same ice sheet sizes at times during 
the previous "80 Kyears? (cf. Andrews, 1973). 

It may be that the seasonal climate part of our 
model is too simplified. Additional non-linearities 
due to land-ocean longitudinal asymmetries, 
realistic atmospheric/oceanic dynamics and struc- 
ture, day-night cycles, cloudiness and pre- 
cipitation variations, etc., could possibly alter the 
model response to the orbital changes to occasion- 
ally give net northern ice sheet budgets of —50 g 
cm”* yr”* (e.g., between 14 and 9 Kyears BP), 
without increasing the amplitude of the intervening 
secondary oscillations of this paper. This possi- 
bility could perhaps best be tested by a higher 
powered GCM, although Hartmann and Short 
0979) and North and Coakley (1979) have shown 
how longitudinal asymmetry can be included 
efficiently in simple climate models. Also Cess and 
Wronka (1979) suggest several new short-term 
feedbacks that could make simple climate model 
response more non-linear. 

Alternatively, our seasonal climate model may 
be basically correct, and there may be other 
long-term processes in the system (apart from ice 
sheet volume inertia) with time scales of several 
Kyears or longer, as proposed by Eriksson (1968), 
Rooth et al. (1978) and others. These processes 
could result in the ice sheet budget depending not 
only on the current ice sheet size but also on its 
past sizes in the previous several Kyears, and so 
could distinguish between the positive and negative 
(—50 g cm”* yr"‘) budgets at the same ice sheet 
size, as discussed above. Specific long-term pro- 
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cesses with this property include time-dependent 
lithospheric depression and/or non-plastic ice sheet 
flow (Emiliani and Geiss, 1957; Tanner, 1965), 
sudden jumps in the profiles of northern hemi- 
spheric and/or Antarctic ice sheets due to basal 
melting or sea incursions (Hollin, 1962; Wilson, 
1964; Hughes, 1977; Thomas and Bentley, 1978), 
and deep ocean temperatures (Newell, 1974; 
Lemke, 1977; Saltzman, 1977). The inclusion of 
two or more long-term processes can produce free 
internal oscillations (e.g., Sergin, 1979), and it may 
be that such an oscillation is important for the 
full glacial-interglacial cycles. 

Two other mechanisms related to the rapid ice 
sheet retreats are a possible layer of ice sheet 
meltwater covering a substantial part of the oceans 
(Adam, 1975; Berger ct al., 1977; Emiliani et al., 
1978), and extensive pro-glacial lakes causing 
calving at the feet of retreating ice sheets (Andrews, 
1973). Processes such as these are not really 
long-term but through them the amount of ablation 
in one year can influence the ice sheet budgets in 
several succeeding years. This could cause a 
climatic flip-flop whereby once large ice sheets start 
to retreat, the oceanic meltwater layer and/or 
proglacial lakes increase and accelerate the retreat 
in succeeding years. 

Before closing, we briefly discuss the relation- 
ship between Saltzman’s (1977) approach and the 
present model. Saltzman shows that the dominant 
balance in the long-term global energy equation 
must be between three terms: global net radiation 
to or from space N, latent heat of fusion associated 
with the changing ice sheet mass F, and thermal 
energy of the global (~deep) oceans W. Newell 
(1974) has described a specific ice age mechanism 
involving these three quantities, and Mason (1976) 
has emphasized the similarity in order of magnitude 
between ice age variations of F and the variation of 


N (at mid and high latitudes) caused directly by th-. 
orbital perturbations. In our model we neglect H 
and our standard method of solution also neglect' 
F, so that N = 0 and global and annual meat 
energy is exactly conserved in eq. (1). Therefore, all 
long-term energy residuals (which are generalh 
much smaller than the seasonal energy terms) art 
neglected in eq. (1), and the long-term changes are 
parameterized by other exploratory equations. Thi^ 
method decreases the computer run-time consider 
ably and has allowed a more extensive exploration 
of diagnostic snow and ice sheet parameteri- 
zations. However, during periods of rapid ice sheet 
retreat the average rate of release of ice sheet latent 
heat of fusion may be ~l/3 that of the insolation 
anomaly at mid and high latitudes due to the orbital 
perturbations (Mason, 1976), and so F should 
perhaps be included in eq. (1). We did include F in 
(1) for one run over the last 100 Kyears and found 
it caused very little difference in the response [a 
nearly constant 2° lat. equatorward shift from the 
solid curve in Fig. 3{b)l. This suggests that the 
“direct” effects of the orbital perturbations on the 
ice sheet budget (via the seasonal climate) are much 
larger than any “adjustments” required to satisfy 
the long-term energy equation, at least for the 
secondary oscillations. 

Relatively simple models such as the present one 
are suitable for experimenting with ice age runs 
incorporating the various mechanisms mentioned 
above. The long-term energy residuals (N, F, and 
IV above) can be explicitly incorporated in the 
model, but are not necessarily important for some 
mechanisms. If any mechanism is found that gives 
realistic glacial-interglacial cycles, the individual 
components and interactions could then be tested 
economically by higher powered GCMs (e.g.. 
Gates, 1976a,b) and ice sheet models (e.g., Jenssen, 
1977). 
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flpH noMoniH npocToR KmiMaTHMecKoil Monejm, 
SBHO onHCbiaaioiueli nenHUKOBbie uiHTbi, Hccnenyexca 
acTpoHOMHHecxafl reopHu jienHHKOBux nepnonoB. 
OflHoypoBeHHaa sonajibHO ocpeAHeHuan cesoHHafi 
3Hepro6anaHCOBaa Monenb peiuaeTcn hhcjichho c 
onpeneJiCHHCM TeMneparypbi na ypoBHC Mopa kbk 
( tiyHKmiH uiupoTbi H Mccflua (noAoOHO Hopcy, 1975). 
CHeroBoil nOKpoB, MeHtnouiHiica c ccsohom (hto 
BAHaer na nnaHerapHoe anbCcAo) BKAioHeH AHarno- 
CTHHCCKH nyrein npocTwx napaMeTpHsauHfl eaceMe- 
caaHbix BbinaAeHHll h TaanHlt cnera. riojiHbie 
roAOBbie HaKoiuieHHa h aOnauHa Ha noBepxHOciH 
neflHHKOBblX UIHTOB H3 KaXCAOli lUHpOTe BblHHCnaiOTCa 
c TCMH ace napaMerpHsauHaMH, mto h ana cnera 
(c TeMnepatypoil, ynHTbiBaioiueli Bwcory lUHTa npH 
BepTHxajibHOM rpaaHCHTe — 6,5°C km"‘). PaccMor- 
peHHC jienuHKOBbix UIHTOB cnenycT YMpTMaHy 
(1976), rae TeacHHe Jibaa annpoKCHMHpyerca xax 
TenenHe HaeanbHO nnacTHMHofi acHAKOCTM, hto 
orpaKHHHBaeT npocJiHJib lUHTa napaSonoii. JleflaHOil 
uiHT ceaepHoro nonyuiapna npocTHpacTca k aKsaxopy 
OT 75°c.m. (mto cooTBexcTByex 6eperoBofl ahhhh 
CeBepHoro jieaoBHXoro oKeana). 


Yapp, C. J. and Epstein, S. 1977. Climatic implications 
of D/H ratios of meteoric water over North America 
(9500-22,000 B.P.) as inferred from ancient wood 
cellulose C-H hydrogen. Earth and Planet. Sci. Lett. 
34, 333-350. 


MoflejibHbie KpHBbie ana nocneaHHx 100 xbic.nex 
paCCHHTaHbl BbIHHCJieHHaMH TOflHHHOii «norOAbl» , 
Kax onHcaHo Bbiuie, ana Kaacawx 2 xbic.aex c hhco- 
aaiiHefl, Bsaxoii c ynexOM peanbHbix BOTMyuieHHH 
napaMCTpOB seMHofl opflHXbi. MsMCHCHHa b pasMepax 
jieaHHKOBoro utHxa ana Kaacaoro uiara b 2 xbic.nex 
aaBHCHT TonbKO OT eaceroaHoro 6anaHca CHcra, 
npoHHxerpHpoBaHHoro no acefi nosepxHOCXH uiHTa. 
B 3THX pacHexax loacHaa rpaHHua neaaHoro uiHxa 
ceaepHoro nonyuiapHa Kone6neTca b npeaenax T 
lUHpOTbl, HTO npaBHAbHO MOaCAMpyeT (Jiasbl H 
npH6nH3HxenbHbie aMnnHxyabi BbicoKoyacxoTHbix 
KOMnoHCHT (43 H 22 Tbic.nex) aaHHbix myOoico- 
BoaHbix ocaaxoB (Xafic h ap., 1976). Oanaxo Moaenb 
HecnocoOna BocnpoHsaecxH ochobhuc neaHHKOBwe 
LUHcnbi (OT 100 ao 120 xwc.nex) sthx aaHHbix. 
OnHcaua HyBCTBHxeJibHOCXb MoaenbHbix pacnexoB 
K H3MeHCHHaM pasAHHHbix ec napaMexpoB, ho hh 
oaHo H3 3THX H3MeHeHMH He ynyHuiaeT cyuiecxBeHHo 
CXOaHMOCTH paCHeXHblX XpHBblX K ynoMaHyxbiM 
aaHHbiM. B saxmoHeHMe mu o6cyacaaeM thhu 
M eXaHH3MOB, KOTOpUe MOTAH 6bl aaXb peaAHCTHHe- 
CKHe neaHHKOBbie umkau. 
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A simple parameterization for ice sheet ablation rate^ 
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ABSTRACT 

A parameterization of monthly mean ice sheet ablation as a function only of surface air 
temperature and insolation is examined. By considering differences in summer climate between 
the present and the last ice age, it is suggested that the parameterization adequately describes 
long-term changes in ablation. Although significant discrepancies are found between the 
parameterization and the seasonal variation of ablation observed on present-day glaciers, the 
seasonal ablation cycle is still predicted accurately enough to maintain the validity of the 
parameterization for long-term net annual variations. 


1. Introduction 

Several recent models of northern hemispheric 
ice sheet fluctuations during the last years 
have used relatively simple parameterizations of 
ablation on the ice sheet surfaces (Weertman, 
1976; Sergin, 1979; Pollard et al., 1980). The 
process of ablation is an important part of the 
climatic control over the size of present-day 
glaciers, and ablation variations in the past could 
have been important in controlling the past ice 
sheet fluctuations. Ablation is used here to mean 
“the reduction in mass of a vertical ice sheet 
column due to the removal of HjO out of the 
column by surface processes”; calving into oceans 
is excluded. In practice ablation mostly involves 
melting and subsequent runoff in surface or basal 
streams during summer, but can also involve 
evaporation or wind-drifting. The full process is 
complex and is described in Paterson (1969, Ch. 4) 
and Sugden and John (1976, Ch. 14). 

In Pollard et al. (1980) the monthly mean 
ablation. A, at any point on the ice sheet surface is 
parameterized in the form 

A — mais\0-, aT + bQ + c\ (1) 

where T is the monthly mean surface air temper- 
ature (corrected for ice sheet elevation above 

' Contribution number 3374 of the Division of 
Geological and Planetary Sciences, California Institute 
of Technology, Pasadena, California 91125, U.S.A. 


sea-level using a constant atmospheric lapse rate), 
Q is the monthly mean insolation at the top of the 
atmosphere, and a, b, and c are constants. For 
non-zero A, (1) is basically a linearized energy- 
balance equation for a melting surface with no heat 
capacity; somewhat more complex parameteri- 
zations involving more climatic variables were first 
developed in the glaciological literature (e.g., 
Sverdrup, 1935). Of course more sophisticated 
snowmelt models are in current use, and a few have 
recently been applied to paleoglaciological 
problems (e.g. Williams, 1979). 

This note examines the adequacy of such a 
simple parameterization as (1) for ice age models. 
In Section 2 past ablation variations for a given 
month (July) are considered, and found to be 
adequately describable in terms of T alone. 
However ablation typically has a large non-linear 
seasonal variation and is negligible in winter, 
making it important to “choose” the right months 
of the year to compare with past eras, i.e., to 
correctly predict the beginning and end of the main 
ablation season. In Section 3 some discrepancies 
are found between present-day seasonal obser- 
vations and a parameterization of the same type as 
(1), but the scatter is only equivalent to a relatively 
small error ('±1 month) in the phase of the 
seasonal ablation cycle. The purpose of including Q 
in (1) and the numerical values for a, b and c are 
discussed in Section 4. 


0040-2826/80/040384-05502.50/0 © 1980 Munksgaard. Copenhagen 


Tellus 32 (1980), 4 



A SIMPLE PARAMETERIZATION FOR ICE SHEET ABLATION RATE 


385 


2. Variations in past eras 

Ablation depends not only on surface air 
temperature and insolation but also on cloud cover, 
wind speed, relative humidity and on physical 
properties of the surface itself. Each of these 
parameters has probably varied systematically over 
synoptic scales in past eras and so could potentially 
be important for ice age ablation parameteri- 
zations. One approach of estimating the relative 
importance of the meteorological parameters (cf 
Coakley, 1977) is shown in Table 1. The ^ values 
in column 2 represent differences in the parameter 
p between the present July and July !8 Kyear BP 
on continents around '-50° N, estimated from the 
GCM results of Gates (1976). These summer 
values are the most relevant to the net annual 
ablation, since by far the most ablation occurs in 
the summer months. The values of the partial 
derivatives BAldp in column 3 are estimated from 
the theoretical “free ablation” graphs of Kraus 
(1975), and represent the regions of these graphs in 
which the bulk of ablation occurs on real glaciers 
(i.e., by melting with air temperature ^0°C); these 
values are consistent with earlier semi-empirical 
ablation formulae (e.g., Sverdrup, 1935). They are 
calculated for a standard glacier surface (with 
albedo 0.5) in “steady-state” ablation, i.e., with no 
sensible heat storages and with all ablated water 
removed immediately. |The partial derivatives for T 


and Q correspond to a and b in (1), and their values 
are considered further in Sections 3 and 4.) 

The magnitudes of the products in column 4 
indicate that by far the largest ablation rate 
changes of past eras have been caused by the 
changes in surface air temperature, T. This sug- 
gests that a parameterization in terms of T alone 
would be adequate for an ice age model; (actually 
the next most important variable, surface wind 
speed, could not realisticaUy be predicted in simple 
one-layer climate models). 

It should be noted that the Ap values in Table 1 
are averages over many days, but the steady-state 
partial derivatives take no account of the pro- 
nounced non-linear daily cycles involved in real 
ablation situations. However in situations where 
melting is occurring for a significant fraction / of 
each day, one would expect that the inclusion of 
these day-night effects would reduce the products 
in column 4 by factors of between '-/ and ~1, and 
would not seriously affect their relative magnitudes. 

Table I also neglects any systematic changes in 
properties of the ice sheets themselves. For instance 
sub-zero internal temperatures can cause melt- 
water to refreeze elsewhere in the ice sheet and not 
run off; the winter cold-wave of seasonal heat in 
near-surface layers can significantly inhibit ablation 
(Muller, 1963). There is little or no direct evidence 
of such systematic changes as ice sheet size varied 
in the past, but these properties do vary widely on 


Table 1 . Effect of past variations Ap of summer climate on ablation rate A 


Parameter p 

Ap (ice age July 
minus present July) 

BAtdp in (g cm"^ month'* 
per units of col. 1) 

Ap (dAidp) in 
(g cm'^ month'*) 

Surface air temperature, T (°C) 

-10 

15 

-150 

Insolation, Q (W m~0 

±15* 

(0.8) X (0.4) = 0.3 2t 

±5 

Total cloudiness (areal fraction) 

+0.2 

Effect of sunlight absorbed: 




(0.8) X (-70) = -56t 

-11 



Effect of net IR: 




(-0.8) X (-80) = 64$ 

+ 13 

Surface wind speed (m s"') 

+2 

15 

+30 

Relative humidity (%) 

+ 10 (at 800 mb level) 

1 

+ 10 


•Representative not of July 18 Kyears BP but of the general magnitude of past variations in mean summer 
insolation at ~50° N due to the orbital perturbations (e.g., Vernekar, 1972). 

+ Factor of (0.8) represents ablation per sunlight absorbed by the glacier surface (Kraus, 1975); factors of (0.4) and 
(—70) represent changes in sunlight absorbed per change in the relevant parameter p as computed by Schneider and 
Dickinson (1976, table I) allowing for multiple reflections between cloud deck and surface. 

t Factor of (—0.8) represents ablation per net infrared radiation lost by the glacier surface to the atmosphere 
(Kraus, 1975); factor of (—80) represents an average sensitivity of this infrared loss to various types of cloud cover as 
determined by glacial field measurements (Wallen, 1948; Lister and Taylor, 1961). 
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present glaciers and ice caps (Paterson, 1969, Ch. 
2). Such properties can be predicted only by 
relatively complex ice sheet models, and are 
effectively assumed constant in most simple ice age 
simulations. 


3. Present-day observations 

In this century considerable glaciological field- 
work has been devoted to relating the ablation on 
glaciers and ice caps to the local meteorological 
conditions (e.g., Paterson, 1969, Ch. 4). Schytt 
(1967) and Loewe (1971) among others have 
compiled sets of present-day glacial data and 
plotted ablation rate. A, against surface air temper- 
ature, T, averaged through most or all of the 
ablation season, i.e., ~May to ^September. (In- 
significant ablation occurs in winter months outside 
of the ablation season.) For interannual variations 
on individual glaciers, the relationship is good and 
shows dAldT -lO (g cm"^ month"') per (°C); 
however there is considerable scatter between 
glaciers in different types of climates (e.g., con- 
tinental versus maritime), and also between 
measurements averaged only over a few weeks or 
less on the same glacier. Much of this scatter in the 
relationship between A and T is probably due to 
seasonal and latitudinal variations of insolation and 
also due to variations in surface albedo, which can 
change seasonally from ^0.8 (fresh snow) to ~0.5 
(melting ice) (e.g., Wallen, 1948). These variations 
would correspond to ablation rate variations of 
~ 100 g cm"^ month in column 4 of Table 1. 

In Fig. 1 we have attempted to improve the 
ablation parameterization by adding a second 
variable “net radiation absorbed”, R, which 
primarily contains the effects of insolation, surface 
albedo and also cloudiness. As far as an accurate 
representation of this data is concerned, the result 
shown is basically negative; the considerable 
scatter in the relationship between A and T has 
been reduced only slightly if at all by the addition 
of the second variable R. The scatter in the figure 
must be due to some combination of observational 
error, seasonal and interglacial variability of ice 
body properties such as seasonal heat storage and 
refreezing, and also due to the exaggerated point- 
to-point variability of surface wind speed and 
relative humidity for measurements averaged only 



Fig. 1. Measurements of present-day ablation vs surface 
air temperature, T, and net solar and infrared radiation 
absorbed, R. Values adjacent to each point are ablation 
rates. 4, in g cm"' month"'. The measurements are 
averaged over a whole number of days during the 
ablation season. The figure contains all the data we found 
where simultaneous determinations of daily mean A, T, 
and R were published. The general level of observational 
error for each variable (not often reported) may be 
<10% of the total range in the figure. The sources for the 
various points in the figure are: • Karsa Glacier 
(Wallen. 1948); + Hoffells. Glacier (Ahlmann and 
Thorarinsson. 1938); * Isaschen. Ice Cap (Sverdrup, 
1935); A Sveanor Snowfield (Sverdrup, 1935); v 
Britannia Glacier (Lister and Taylor, 1961); x Barnes 
Ice Cap (Ward and Orvig, 1953); □ Kessel. Snowfield 
(Ambach and Hoinkes, 1963). The dashed lines of 
constant ablation rate show a “fit" that corresponds to 
the parameterization in the text. 

over a few days. Unfortunately even if these effects 
were resolvable from the available data, a realistic 
treatment of these variables (and also of season- 
ally varying surface albedo contained in the 
variable R) is beyond the scope of most simple ice 
age models. 

The level of scatter from the linear fit in Fig. 1 is 
<±50 g cm"' month"'. This is equivalent to 
--±5°C in air temperature, which is much less 
than the full seasonal variation at any one location 
and is comparable to the change through one 
month in spring or autumn. Therefore the begin- 
ning and end of the main ablation season could be 
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predicted to within month by a parameteri- 
zation in terms of T alone, with a non-linear cutoff 
below T^—5 °C to represent negligible ablation in 
winter. 

4. Parameterization 

Although the radiative variable R has not 
eliminated the scatter in Fig. 1, it still seems 
preferable to include the seasonal variation of Q in 
an ablation parameterization for the following 
reason: with all other variables held constant, the 
full winter-summer variation of Q would certainly 
affect ablation as much as a variation in T of 
~10°C, and so the inclusion of Q should slightly 
improve the predicted phase of the seasonal cycle. 

Numerical values for (1) are chosen below, but 
one should note that these values are constrained 
only within wide limits by the data in this paper, 
and the main point is the basic form of (1). As in 
Table 1, we estimate 8A/dQ ~0.3 (g cm“* 
month” ') per (W m”*). (This steady-state 
derivative should be affected only slightly by 
day-night effects since the daily cycle of insolation 
is roughly correlated with that of ablation.) Using 
the value of dA/dT from present-day interannual 
variations (Section 3), this yields the following 
equation for monthly mean quantities: 

A (g cm”^ month”') s; max[0; lOF (°C) 

, -t- 0.3(2 (Wm-^)- 50] (2) 


We have estimated the constant “50” in (2) by 
comparison with the data in Fig. 1, using R (W 
m“^) = 0.4g — 80. [The factor 0.4 allows for 
multiple reflections off a ~0.3 cloud cover 
(Schneider and Dickinson, 1976), and 80 W m”^ 
represents net infrared loss from a melting surface 
with this cloud cover (Lister and Taylor, 1961).] 
The corresponding linear fit is shown in the figure. 
This constant “50” in (2) is also consistent with the 
general observation that significant ablation occurs 
only where the monthly mean temperature rises 
much above - — 7°C in summer (e.g., Orvig, 1954; 
Bull and Carnein, 1968). The non-linear cutoff in 
(2), basically due to the transition from melting to 
evaporation, corresponds to the non-linear temper- 
ature dependence (roughly to T^) suggested by 
Ahlmann (1948). 

In summary, one should expect significant 
discrepancies between eq. (2) and present observed 
seasonal cycles of ablation. Much of this dis- 
crepancy could be due to seasonal variations of 
surface albedo, but Fig. 1 shows that other factors 
(e.g., seasonal heat storage) are significant. How- 
ever, any such seasonal discrepancy should not be 
serious for the ice age problem as long as the 
summer months during which most ablation occurs 
are predicted reasonably well. This is because 
Table 1 shows that for these months the past 
variations of ablation rate have been controlled by 
past variations of surface air temperature, con- 
sistent with eq. (2). 
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EMPIRICAL RELATIONS BETWEEN RADIATIVE FLDl, 
CLOUDINESS. AND OTHER METEOROLOGICAL PARAMETQtS 


D.D. Wenkert and A.P. Ingertoll 
California Inatitote of Technology 
Paaadena, California 


1. INTRODUCTION 

We are analyzing Earth Radiation Bndget 
(ERB) data obtained by the scanning, narrow 
f ield-of-viev channels of the ERB instrament on 
Nimbus 7 from the whole Earth for six days in 
June 1979. Using meteorological information froai 
the FGGE Level Illb (ECHWF) data set and cloudi- 
ness information from the THIR instrument onboard 
Nimbus 7. we calculate statistical relations 
between the angularly resolved radiances measured 
by ERB and various meteorological and cloud 
parameters. These relations are then combined to 
determine relations between angularly integrated 
visible and infrared upwelling flux (irradiance) 
and meteorological parameters (including 
cloudiness) . 

We ate using the relationships we 
obtain between radiative flux and other atmos- 
pheric variables in three ways. The effect of 
varying cloudiness on the net flux at the top of 
the atmosphere is being determined. Various dif- 
ferent radiation parameterizatlons. appropriate 
for simple climate models, are being developed. 
Finally, relationships between radiation and 
other atmospheric variables are being defined, 
against which the internal statistics of general 
circulation models may be tested. In this paper, 
we report on some preliminary results oriented 
primarily toward solving the first problem. 

2. OBSERVATIONS 

The radiation data consist of angularly 
resolved, broad band measurements of infrared and 
visible radiation obtained by the four scanning 
telescopes of the ERB instrument on the Nimbus 7 
satellite, which travels in a sun-synchronous 
noon-midnight orbit. At any given time, the 
scanners were capable of sampling most of the 
Earth's visible disc with instantaneous fields of 
view of 0.2S degrees x 5.12 degrees (Jscobowitz 
et al.< 1978). The data are archived on Sob 
Target Radiance Tapes (STRT's, tee Stowe and 
Fromm, 1983), where they are organized by target 
area (TA, a roughly square region on Earth about 
500 hm on a side; there are 2070 of them), sob 
target area (there are nine of these in each TA), 
time of day, and viewing geometry. Also included 
on the STRT are information on the type of 
surface found in each TA and STA, and cloudiness 


information for the STA'a (when available) 
obtained by the THIR instrument on Nimbus 7 (Chen 
et al.. 1980). 

The meteorological information we use 
comet from the 6-honrly global synoptic analyses, 
prodnced at a Level Illb data set by the European 
Center for Medlnm-range Weather Forecasting 
(ECMWF) for the second Special Observing Period 
of FGGE. This data set hat space and time reso- 
Intiona that match quite well those of the ERB 
data. 

3. STATISTICAL ANALYSIS 

In general, the upwelling flux, F, from 
a region it related to the upwelling radiance 
field. R. by 

P(po) ) i R(0.d;Mo) cos 6d 0 sin 6 dd. 

2n 

where p, « cosine of tolar zenith angle, 6 = 
zenith angle of observer, i = relative azimuth 
from snn to observe. We have divided the upward 
hemisphere into 32 discrete angular bins for the 
visible radiances and 14 angular bins for the 
infrared. Thus the relation between flux and 
radiance becomes 

F(|*c) = Ri(|io) Pi Ci . 

i 

where ■= cosine of zenith angle in ith bin, 

0^ = size of ith angular bin in steradians. It 
is simple to show, given any variable x. that 

"iF “ Z ®xRi “i “i 
1 

where = covariance of x and F (flux), = 

covariance of x and R^ (radiance in 1th bi^. 
Similarly. 

®F^ “ Z Z ®RiRj t*i Oi l*j Oj 
1 j 

where Op = variance of flux cr^^pj = covariance 
of radiance in ith bln and radiance in jth bin. 
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BESULTS 


It is rot necessary to neasore the entire 
upwelling radiation field each tine to deternine 
statistical relationships between upwelling flux 
and other parameters. All that is necesaary is 
that covariancea of the parameter in qneation 
with the radiances of all the angular blna be 
obtained and that covariances of the radiance in 
each angular bin with the radlancea in each of 
the other bins be known. If at least a few times 
during the period under investigation the radi- 
ance in each bin has been measured at the same 
time (i.e., within the same time averaging inter- 
val) as the radiance in any given other bin. and 
if the same is true for the radiance in each bin 
and each of the meteorological (and cloudiness) 
parameters in which we are Interested, then we 
have sufficient information to deduce statistical 
relations between the upwelling radiative flux 
and the parameters of interest. 

For the investigations discussed here 
we are interested in the dinrnally averaged 
radiative flux rather than the flux at one 
specific solar zenith angle. Since the ERB 
instrument takes data only near local noon and 
local midnight, this can be a problem for visible 
radiance measurements. If we let x be the 
diurnal average of x. and let t be the time of 
day then 

F(pa(t)) = 5lRi(i*a(t)) . 

and F “ 

i 

We can define a bidirectional reflectance such 
that 

Ej(t) = Mod) Aj^Ct) Fo/n . 

where Fo tolar flux (solar constant), 

A^(t) '= A£(Mo(t),s) ■= bidirectional reflectance 
for a cloud and surface type t, for a cosine of 
solar zenith angle Mo> d the ith outgoing 
angular bin. Now we can use 

R"i = MoAj Fo/« = (m^A^/mo ( t) A^(t)) R£(t) 

If we know the way in which the bidirectional 
reflectance varies with solar zenith angle then 
it is no problem to determine the diurnally 
averaged visible radiance from an observed 
instantaneous visible radiance. For the results 
reported here, we have assumed that in all cases 
A^ was independent of po* Thus we used 

= (pi/po(t)) Rj^(t) 

Note that p^ it a known function of latitude and 
time of year and Po(t) is listed in the SIRT 
records for all observations. The averaging 
interval we use for our data is 24 hours and the 
averaging region is a full Target Area (c. 500 km 
square). However, we do not average in any 
visible radiances measured for regions in 
darkness. It is hoped that the scheme outlined 
above provides good diurnal averages for visible 
radiances and that averaging both the near noon 
and near midnight values of infrared radiances 
provides good diurnal averages for those. Using 
standard linear multiple regression, the vari- 
snces and covariances of radiative fluxes and 
other parameters yield empirical equations 
relating these variables to each other. 


4 . 

In Table 1, we have liated the coeffi- 
cients for single parameter regressions of the 
different (dinrnally averaged) radiative fluxes 
against 1000 mbar temperature and against 
fractional cloudiness (as determined . by THIR). 

The net flux (Fnet) is simply the average inci- 
dent tolar flux minus the sum of the outgoing 
visible and Infrared fluxes. The regressions 
were determined from the data set for one day (12 
June 1979) for the varioua regions of the Earth 
listed. The single parameter regression coeffi- 
cients are equivalent to total derivatives of the 
flux with respect to the appropriate parameter 
(either TIOOO or Xcld). Thus the variation in 
radiation doe to cloudiness or surface air tem- 
perature, listed in this table reflect the actual 
distributions of cloudiness or surface air tem- 
perature on 12 June 1979. To determine the 
effect, for example, of varying cloudiness while 
holding latitude and surface air temperature 
constant, a mnltlpsrameter regression should be 
performed. In order to get an idea of the magni- 
tude of the different variations in radiation doe 
to a change in cloudiness or temperature, we have 
listed, in Table 2, the change in flux due to a 
positive change in TIOOO or Xcld of one standard 
deviation. These are compared to the standard 
deviations of the fluxes. It is clear from these 
tables, that on 12 June 1979 higher cloudiness 
was associated with effects on outgoing terres- 
trial radiation and reflected solar radiation 
that nearly cancel when averaged over the entire 
Earth. This near cancellation is also seen in 
the data sets for most of the other five days 
investigated, and it shows up in mnltlparameter 
regressions Involving cloudiness. Therefore we 
believe this effect is real. The breakdown of 
the cloudiness effect by latitude band it 
enlightening. The effect of clouds on emission 
dominates in those regions receiving little inso- 
lation (i.e., the winter hemisphere) and the 
albedo effect dominates in regions receiving more 
sunlight (i.e. the summer hemisphere). The net 
result over the globe is near cancellation. 
Increased surface air temperature is associated 
with decreased net flux over most of the globe, 
presumably because of increased thermal emission 
with increasing temperature. Interestingly, 
however, this relationship breaks down in some 
regions, specifically the tropics. Perhaps in 
this region. Increasing surface temperatures lead 
to increasing atmospheric 

humidity and penetrative convective cloudiness, 
both leading to a higher atmospheric opacity. 


5. (X)NCLUDING REMARKS 

The Nimbus 7 ERB data set can be used 
along with correlative data to investigate the 
effects of varying cloudiness and meteorology on 
the radation budget of the atmosphere, without 
using any a priori assumptions or angular models. 
However, to properly correct for the probable 
systematic underestimation of dinrnally averaged 
upwelling visible flux, due to oollection of 
visible ERB data at the lowest tolar zenith 
angles (corresponding to noontime) for each 
target, tome assumptions about angular models 
most be made. We are presently developing 
various diurnal corrections which involve a mini- 
mum of assumptions. We are also investigating 
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the effecte of differeot averaging intervala (of 
area and tine) and different divitions of the 
upward henitphere Into angnlar bina on onr atatl- 
atica. Ve hope to thereby nahe nore robnat onr 
findings on the relationships between clondlneaa, 
neteorology. and atnospherlc radiation, inclndlng 
the apparently neutral effect of clouds on the 
Earth's radiation budget. 
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Table 1 


Coefficients for single psrsiieter regressions 


Region 

Parameter 

Fvis 

Fir 

Fnet 

90. OS-67. 5S 

TIOOO 

O.OOOE-sOO 

3.205E-K>0 

-3.205E+00 


Xcld 

O.OOOE-^OO 

-5.881E-01 

5.881E-01 

67.5S-45.0S 

TIOOO 

2.394E-K)0 

1.455E+00 

-3.849E-K)0 


Xcld 

1.552E+01 

-6.315E+01 

4.763E+01 

45. OS-22. 5S 

TIOOO 

5.914E-01 

2.843E+00 

-3.435E+00 


Xcld 

4.608E-^01 

-8.768E+01 

4.160E-)-01 

22.5S-0.0 

TIOOO 

l.dllE-fOO 

-4.457E+00 

3.045E-t-00 


Xcld 

6.335E-»^01 

-6.589E-t-01 

2.541E+00 

0.0-22.5N 

TIOOO 

-9.079E-01 

-2.590E+00 

3.498E-I-00 


Xcld 

8.097E+01 

-1.321E+02 

5.111E+01 

22.5N-45.0N 

TIOOO 

-2.781E+00 

3.571E+00 

-7.902E-01 


Xcld 

7.182E+01 

-1.024E+02 

3.060E-S01 

45.0N-67.5N 

TIOOO 

-7.308E+00 

2.204E-M)0 

5.104E-^00 


Xcld 

1.870E-S02 

-7.299E+01 

-1.140E+02 

67.5N-90.0N 

nooo 

-4.846E+00 

5.602E-01 

4.286E-S00 


Xcld 

2.026E402 

-6.178E+01 

-1.409E+02 

Whole Earth 

TIOOO 

6.241E-01 

2.552E+00 

-3.176E+00 


Xcld 

8.991E+01 

-8.402E+01 

-5.884E+00 


Table 2. 

Effects on radiative fluxes of ebangingtbe regression 
parameter by standard deviation froai the mean 
(standard deviations of the fluxes included for coaiparlton) 


Region 

Parameter 

Fvis 

Fir 

Fnet 

90. OS-67. 5S 

TIOOO 

O.OOOE+00 

1.570E+01 

-1.570E+01 


Xcld 

O.OOOE+00 

-1.907E-01 

1.907E-01 

Standard 

Deviations: 

O.OOOE+00 

2.503E-)-01 


67.5S-45.0S 

TIOOO 

1.320E+01 

8.026E+00 

-2.123E+01 


Xcld 

2.611E+00 

-1.062E+01 

8.010E+00 

Standard 

Deviations: 

1.476E+01 

2.138E-M>1 


45. OS-22. 5S 

TIOOO 

2.614E-^00 

1.257E+01 

-1.518E+01 


Xcld 

1.201E+01 

-2.285E+01 

1.084E-M)1 

Standard 

Deviations : 

1.853E+01 

2.775E+01 


22.5S-0.0 

TIOOO 

3 .250E+00 

-1.026E+01 

7.012E+00 


Xcld 

1.777E+01 

-1.849E+01 

7.128E-01 

Standard 

Deviations : 

2.591E+01 

2.513E+01 


0.0-22.5N 

TIOOO 

-1.790E+00 

-5.107E+00 

6.897E-)-00 


Xcld 

1.676E+01 

-2.733E+01 

1.058E-H)1 

Standard 

Deviations: 

4.186E-*-01 

4.365E+01 


22.5N-45.0N 

TIOOO 

-1.050E+01 

1.348E-^01 

-2.983E+00 


Xcld 

2.152E+01 

-3.069E+01 

9.169E-f00 

Standard 

Deviations: 

3.793E+01 

3.833E-*-01 


45.0N-67.5N 

TIOOO 

-3.806E+01 

1.148E+01 

2.658Em 


Xcld 

4.278E+01 

-1.670E+01 

-2.608E+01 

Standard 

Deviations : 

5.863E+01 

2.290E-K01 


67.5N-90.0N 

TIOOO 

-2.277E+01 

2.633E+00 

2.014E-^01 


Xcld 

2.266E+01 

-6.910E+00 

-1.575E+01 

Standard 

Deviations : 

4.099E+01 

1.071E+01 


Whole Earth 

TIOOO 

6.793E+00 

2.778E+01 

-3.457E+01 


Xcld 

2.416E-S01 

-2.258E-K)1 

-1.581E+00 

Standard 

Deviations: 

6.732E+01 

4.324E->^01 
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